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Abstract 

The UNEDF project was a large-scale collaborative effort that applied high-performance computing to the nuclear 
quantum many-body problem. UNEDF demonstrated that close associations among nuclear physicists, mathemati- 
cians, and computer scientists can lead to novel physics outcomes built on algorithmic innovations and computational 



^^f^ developments. This review showcases a wide range of UNEDF science results to illustrate this interplay. 
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Understanding the properties of atomic nuclei is cru- 
cial for a complete nuclear theory, for element formation, 
for properties of stars, and for present and future energy 
and defense apphcations. From 2006 to 2012, the UNEDF 
(Universal Nuclear Energy Density Functional) collabo- 
ration carried out a comprehensive study of the nuclear 
many-body problem using advanced numerical algorithms 
and extensive computational resources, with a view toward 
scaling to petaflop supercomputing platforms and beyond. 

The UNEDF project was carried out as part of the 
SciDAC (Scientific Discovery through Advanced Comput- 
ing) program led by Advanced Scientific Computing Re- 
search (ASCR), part of the Office of Science in the U.S. 
Department of Energy (DOE). The SciDAC program was 



started in 2001 as a way to couple the applied mathemat- 
ics and computer science research sponsored by ASCR to 
applied computational science application projects tradi- 
tionally supported by other offices in DOE. UNEDF was 
funded jointly by ASCR, the Nuclear Physics program of 
the Office of Science, and the National Nuclear Security 
Administration. Over 50 physicists, applied mathemati- 
cians, and computer scientists from 9 universities and 7 
national laboratories in the United States, as well as many 
international collaborators, participated in UNEDF. 

This review describes science outcomes in nuclear 
many-body physics, with an emphasis on computational 
and algorithmic developments, that have resulted from 
the successful collaborations within UNEDF among math- 
ematicians and computer scientists on one side and nuclear 
physicists on the other. Such collaborations "across the di- 



vide" were newly formed at the early stage of the project 
and became its unique feature, with high-performance 
computing serving as a catalyst for new interactions. 
The results described in this paper could not have been 
achieved without such couplings. 

1.1. UNEDF science 

The long-term vision initiated with UNEDF is to ar- 
rive at a comprehensive, quantitative, and unified descrip- 
tion of nuclei and their reactions that is grounded in the 
fundamental interactions between the constituent nucle- 
ons [U 12] . The goal is to replace phenomenological models 
of nuclear structure and reactions with a well-founded mi- 
croscopic theory that delivers maximum predictive power 
with well- quantified uncertainties. Specifically, the mis- 
sion of UNEDF was threefold: 

1. Find an optimized energy density functional (EDF) 
using all our knowledge of the nucleonic Hamiltonian 
and basic nuclear properties. 

2. Validate the functional using the relevant nuclear 
data. 

3. Apply the validated theory to properties of interest 
that cannot be measured. 

The main physics areas of UNEDF, defined at the be- 
ginning of the project [1 , were ab initio structure, ab ini- 
tio functionals, density functional theory (DFT) applica- 
tions, DFT extensions, and reactions. Few connections 
between these areas existed at that time. As UNEDF ma- 
tured, however, coherence grew within the effort. Indeed, 
the project created and facilitated an increasing interplay 
among the major areas where none had existed previously. 
Each of the main physics areas now includes ongoing col- 
laborations that cross over into other areas. These in- 
terconnections are highlighted in the summary diagram 
of the UNEDF strategy shown in Fig. [l] In addition to 
physics links, numerous computer science/applied mathe- 
matics (CS/AM) interconnections were established within 
UNEDF as computational and mathematical tools devel- 
oped in one area of UNEDF were used in other parts of the 
project. These tools, motivated by nuclear needs, are now 
available for other areas of science. Access to leadership- 
class computing resources and large-scale compute time 
allocations were critical for the scientific investigations. 

At the intersection of the ab initio techniques and DFT 
techniques are comparisons of observables among the var- 
ious approaches, particularly through constraints on den- 
sity. Such calculations have not been performed before 
and require significant computational capability and an in- 
creasing sophistication of data manipulation. Research on 
the nuclear problem would be incomplete without a serious 
effort to understand the nuclear interactions involved and 
their connection to DFT. Therefore, the UNEDF project 
also included elements that required less computational 
capability but are integral to the project, such as the de- 
velopment of nuclear forces using renormalization group 
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Figure 1: UNEDF project scope. Major science areas are indicated 
by boxes; interconnections between areas are marked by arrows. The 
green boxes indicate connections to experimental observations. 



approaches. Another example is research on nuclear reac- 
tion properties that requires both the use and development 
of algorithms for the largest computers and more conven- 
tional computing needed for algorithmic breakthroughs. 

Another new aspect of the nuclear theory effort driven 
by this project is a greatly enhanced degree of quality con- 
trol. Integral to UNEDF was the verification of methods 
and codes, the estimation of uncertainties, and other out- 
put assessments. Methods used for verification and val- 
idation included the crosschecking of different theoreti- 
cal methods and codes, the use of multiple DFT solvers 
with benchmarking, and benchmarking of different ab ini- 
tio methods using the same Hamiltonian. A new way to es- 
timate theory error bars was to use multiple Hamiltonians 
with different energy/momentum cutoffs and then analyze 
the cutoff dependence of calculated observables. The UN- 
EDF assessment component necessitated the development 
and application of statistical tools to deliver uncertainty 
quantification and error analysis for theoretical studies as 
well as to assess the significance of new experimental data. 
Such technologies are essential as new theories and compu- 
tational tools are applied to entirely new nuclear systems 
and to conditions that are not accessible to experiment. 



1.2. Collaborative effort 

The successes of the UNEDF project were built upon 
certain best practices, some implemented originally and 
some learned by experience, in organizing and implement- 
ing the scientific effort. In order to foster the close align- 
ment of the necessary applied mathematics and computer 
science research with the necessary physics research, mul- 
tiple direct partnerships were formed consisting of com- 
puter scientists and applied mathematicians linked with 
specific physicists to remove algorithmic and/or computa- 
tional barriers to progress. The five-year lifetime of the 
project provided time for these collaborations to become 
deep, and they have continued into follow-up projects. 

All these partnerships have success stories to tell, from 
greatly improved load balancing on leadership-class ma- 
chines, to new DFT solver technologies, to dramatically 
improved algorithms for optimization of functionals, to 
eigenvalues and eigenfunctions of extremely large matri- 
ces, and more. 

The SciDAC program aims at transformative science, 
and this goal has been fulfilled by the new capabilities 
stemming from UNEDF. But the outcomes reach beyond 
the many compelling nuclear physics calculations. UN- 
EDF has changed for the better the way that low-energy 
nuclear theory is carried out, analogous to the shift in 
experimental programs, moving from many small groups 
working independently to large-scale collaborative efforts. 

2. Science 

The territory of UNEDF science is the chart of the nu- 
clides in the (A^, Z)-plane shown in Fig. |2] On this chart, 
stable nuclei are represented by black squares, while the 
yellow squares indicate unstable nuclei that have been seen 
in the laboratory. The sizable green area marked "terra 
incognita" is populated by unstable isotopes yet to be ex- 
plored. Above the table of nuclides are shown three broad 
classes of theoretical methods, which are also used in other 
fields dealing with strongly interacting many-body sys- 
tems, such as quantum chemistry and condensed matter 
physics. Light nuclei and their reactions can be computed 
by using ab initio techniques (quantum Monte Carlo, no- 
core shell model) described in Sec. 2.1 Medium-mass nu- 



clei can be treated by configuration interaction (CI) tech- 
niques (Sec. |2.2[ ). The bulk of the nuclides are covered by 
the nuclear DFT described in Sec. |2.3[ which provides the 
theoretical underpinning and computational framework for 
building a nuclear EDF. Time-dependent phenomena in- 
volving complex nuclei, including nuclear reactions, can 
be described by means of approaches going beyond static 



DFT (Sec. 2.4). By enhancing and exploiting the overlaps 
with ab initio and CI approaches, the goal is to construct 
and validate a nuclear EDF informed by microscopic in- 
teractions as well as experimental data. 
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Figure 2: Theoretical approaches for solving the nuclear quantum 
many-body problem used by UNEDF. The lightest nuclei can be 
computed by using ab initio methods based on the bare internucleon 
interactions (red). Medium-mass nuclei can be treated by configu- 
ration interaction techniques (green). For heavy nuclei, the density 
functional theory based on the optimized energy density functional 
is the tool of choice. (From J^.) 



2.1. Ab initio methods and benchmarking 

Ab initio methods solve few- and many-body prob- 
lems by using realistic two- and three-nucleon interactions 
and obtain the structure and dynamic properties of nuclei. 
The nuclear interaction depends on the spatial, spin, and 
isospin coordinates of the nucleons. Consequently, calcu- 
lations are much more computationally demanding than 
typical quantum problems. Items of interest include nu- 
clear spectra, charge and magnetic ground-state and tran- 
sition densities, electron and neutrino scattering, and low- 
energy reactions. The main goals are to reproduce known 
nuclear properties and predict properties that are difficult 
or impossible to measure. 

Several ab initio methods have been developed for 
studying light nuclei; all have analogues in the study 
of condensed matter and electronic systems. Quantum 
Monte Carlo (QMC) methods, including Green's func- 
tion Monte Carlo (GFMC), use Monte Carlo evaluations 
of path integrals, explicitly summing over the spin states 
and isospin states of the system. The most recent GFMC 
calculations have concentrated on the ^^C nucleus, a fas- 
cinating system with a low-lying excited 0+ state, the 
Hoyle state, very near the threshold of three- alpha par- 
ticles. QMC methods have also been used to calculate the 
properties of neutron matter and neutrons in inhomoge- 
neous potentials. 

No-core shell model (NCSM) methods, including the 
large-scale many-fermion dynamics nuclear (mfdu) code, 
expand the interacting states in products of single-particle 
states and project the low-lying states through large-scale 
matrix operations. MFDn calculations have been used, for 
example, to explain the long lifetime of the ^^C nucleus 



used in carbon dating. A combination of no-core shell 
model techniques with the resonating group method is cur- 
rently used to calculate important low-energy nuclear re- 
actions. 

The coupled-cluster method is an ideal microscopic ap- 
proach to describe nuclei with closed (sub) shells and their 
neighbors. It exhibits a low computational cost (scales 
polynomially with system size) while capturing the dom- 
inant parts of correlations in the wave function. This 
method has been employed to describe and predict the 
structure and reactions of neutron-rich oxygen and cal- 
cium isotopes. 
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Figure 3: Weak scaling of AGFMC with ADLB in terms of MPI ranks. 
There are 8 ranks per BG/Q node; each rank is using 6 OpenMP 
threads. Note the compressed vertical scale. 



2.1.1. GFMC 

Green's function Monte Carlo calculations start with 
an initial trial state ^t and obtain expectation values in 
the exact eigenfunction ^o of the Hamiltonian. These 
calculations are done by evolution in imaginary time r: 
^0 = exp[— J^rJ^T for sufficiently large r. The evolution 
is done in many small steps of r, each step being a nested 
3A-dimensional integral. GFMC was introduced in light 
nuclei ^ [4] to include the strong correlations induced by 
the nuclear interaction. This method has been used to 
calculate the spectra of light nuclei up to ^^C [H [5], as 
well as form factors, electron scattering, and low-energy 
reactions [6 . 

Calculations of ^^C require the largest-scale computers 
available, using a combination of efficient load-balancing 
for the Monte Carlo and large-scale linear algebra for the 
spin-isospin degrees of freedom. The calculations of ^^C 
required the development of the Asynchronous Dynamic 
Load Balancing (adlb) library to efficiently perform the 
load balancing on more than 100,000 cores [F. 

A program, AGFMC, has been developed over the past 
15 years to carry out these calculations [ll [3 [9]. It is a 
large (80,000 lines) Fortran code that originally used MPI 
to manage parallelism. At the beginning of this project, 
the AGFMC code was scaling well up to around 2,000 pro- 
cesses and performing satisfactorily on IBM's Blue Gene/L 
computer. At that time it was becoming apparent that if 
the code were to be able to take advantage of new, petas- 
cale machines expected to come on line during the five-year 
project to investigate larger nuclei, a significant increase 
in the degree of parallelism would need to be incorporated 
into its main algorithms. The greater degree of parallelism 
(from thousands to tens of thousands of processes) would 
give rise to load-balancing problems that would strain the 
then-used approach. 

One of the goals of UNEDF was to construct a soft- 
ware library, intrinsically general-purpose but with fea- 
tures driven by the requirements of AGFMC, to attack the 
load-balancing problem. The purposes of the library were 
to supply a programming interface that would enable rel- 
atively straightforward migration of the existing AGFMC 
code to the new load-balancing library and to scale the 
entire system to much larger degrees of parallelism. 



The result is the adlb library [5 . adlb generalizes the 
classical manager-worker parallel programming model by 
allowing application processes (workers) to put arbitrary 
independent work units into a shared pool and get them 
out to complete them, notifying other processes when they 
have done so. Work units are assigned types and priorities 
by the workers and retrieved according to these proper- 
ties, allowing complex algorithms to be implemented, de- 
spite the simple nature of the parallel programming model. 
Scalability is achieved by dedicating a small percentage 
(but still potentially a large number) of the job's pro- 
cesses to maintaining this work pool and responding to 
put and get requests. These "server" processes execute 
independently from the application processes, thus allow- 
ing asynchronous load balancing of process load, memory 
consumption for the work pool, and message traffic. 

This scheme has worked well. Most of the MPI pro- 
gramming in the original AGFMC code has been absorbed 
into the adlb library, yet the overall code structure has 
been maintained. Scalability has been extended to more 
than 32,000 processes on BG/P and more than 260,000 
processes on BG/Q (see Fig.pl, enabling scientific results 
unattainable before this project was undertaken. 

The ^^C nucleus is particularly intriguing because it 
has a low-lying 0+ excited state (the "Hoyle" state) very 
near the energy of the breakup into three alpha particles. 
This state is essential for the nucleosynthesis of carbon in 
stars through the triple-alpha process. For ^^C the ^t 
are linear combinations of shell-model and alpha-cluster 
states. Figure [4] shows the convergence of the calcula- 
tions of the ground and Hoyle states in the AGFMC calcula- 
tions. Two different sets of initial states are propagated to 
r ~ l.OMeV" ; they yield consistent results. The ground- 
state energy is well reproduced, and the Hoyle state ex- 
citation energy is approximately reproduced (see [TOl [11] 
for complementary calculations of the Hoyle state). The 
ground-state form factor of ^^C is also reproduced by these 
calculations. 

Other recent applications of AGFMC include pair mo- 
mentum distributions [12], electromagnetic transitions 
[13 , and the studies of trapped neutrons ("drops") de- 
scribed in Sec. 12.3.41 
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2.1.2. NCSM and MFDn 

The measured lifetime of ^^C, 5730±30 years, is a valu- 
able chronometer for many practical applications ranging 
from archeology to physiology. It is anomalously long com- 
pared with lifetimes of other light nuclei undergoing the 
same decay process, allowed Gamow- Teller (GT) beta de- 
cay. This lifetime poses a major challenge to theory be- 
cause traditional realistic nucleon-nucleon (NN) interac- 
tions alone appear insufficient to produce the effect [14 . 
Since the transition operator, in leading approximation, 
depends on the nucleon spin and charge but not the spatial 
coordinates, this decay provides a precision tool to inspect 
selected features of the initial and final nuclear states. To 
convincingly explain this strongly inhibited transition, we 
need a microscopic description that introduces all physi- 
cally relevant 14-nucleon configurations in the initial and 
final states and a realistic Hamiltonian. 

Since the nuclear strong interaction governs the config- 
uration mixing, the Hamiltonian matrix eigenvalue prob- 
lem is a very large, sparse matrix in the configuration space 
of 14 nucleons. We address this computational challenge 
with the MFDn code ^ Ull UHl |19J . Aided by a collabora- 
tion with applied mathematicians on scalable eigensolvers 
and computational resources on leadership-class machines, 
we are able to solve this beta decay problem with sufficient 
accuracy to resolve the puzzle: the decay is inhibited by 
the role of 3-nucleon forces (3NFs) as shown in Fig. [s] (see 
pQ] for complementary calculations). 

We obtained our results on the Jaguar supercomputer 
(see Sec. [4]) using up to 35,778 hex-core processors (214,668 
cores) and up to 6 hours of elapsed time for each set of 
low- lying eigenvalues and eigenvectors. The number of 
nonvanishing matrix elements exceeded the total memory 
available and required matrix element recomputation "on 
the fly" for the iterative diagonalization process employing 
the Lanczos algorithm. 

These calculations and many other achievements [21 
were made possible by dramatic improvements to MFDn 
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described by a chiral effective field theory interaction (adopted from 
15 ). The top panel displays the contributions with (two right bars 
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at A^max = 8. Contributions are summed within each shell to yield 
a total for that shell. The bottom panel displays the running sum of 
the GT contributions over the shells. Note the order-of-magnitude 
suppression of the 0j9-shell contributions arising from the 3NFs. 



capabilities during the UNEDF project [22]. The current 
scaling performance of MFDn is demonstrated in Fig. [6] 
Other recent applications of MFDn include the prediction 
(before experimental confirmation) of the spectroscopy of 
proton-unstable ^^F [23] and studies of trapped neutrons 
( "drops" ) with a variety of interactions and other ab initio 
computational methods [24] . 

2.1.3. NCSM and the resonating group method 

Weakly bound nuclei, or even unbound exotic nuclei, 
cannot be understood by using only bound-state tech- 
niques. Our ab initio many-body approach, no-core shell 
model with continuum (NCSMC), focuses on a unified 
description of both bound and unbound states. With 
such an approach, we can simultaneously investigate struc- 
ture of nuclei and their reactions. The method com- 
bines square-integrable harmonic-oscillator basis (i.e., via 
the NCSM [21 ) accounting for the short- and medium- 
range many-nucleon correlations with a continuous basis 
(i.e., via the NCSM with the resonating group method 
(NCSM/RGM) [25, 26 ) accounting for long-range corre- 
lations between clusters of nucleons. With this technique, 
we can predict the ground- and excited-state energies of 
light nuclei (p-shell, A<1Q) as well as their electromagnetic 
moments and transitions, including weak transitions. Fur- 
thermore, we can investigate properties of resonances and 
calculate characteristics of binary nuclear reactions (e.g.. 
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Figure 6: Strong scaling for MFDn: speedup for 500 Lanczos itera- 
tions (the most time-consuming phase of the code). Two problems 
are shown with their dimension (D) and number of nonzero matrix 
elements (NNZ) in the legend. The smaller is ^Li (D=6.2 million, 
NNZ=118 billion), and the larger is ^^B (D=160 million, NNZ=5.2 
trillion). The smaller problem needs at least 1 TB in order to store 
all nonzero matrix elements in core and needs, therefore, at least 728 
cores to fit the problem in core. The larger problem needs at least 42 
TB, and we used between 30,624 and 261,120 cores for that problem. 



Figure 7: Experimental results for S-factor of ^iie{d,p)'^¥Le reaction 
from beam-target measurements. The full line represents the ab ini- 
tio calculation. No low-energy enhancement is present in the theoret- 
ical results, contrary to the laboratory beam-target data represented 
by symbols; see 28 for details. 



cross sections, analyzing powers). 

Recent applications of our ab initio techniques include 
an investigation of the unbound ^He [27 , calculations of 
^H((i,n)'^He and ^'Re{d,p)^Re fusion [28 (see Fig.[7|), and 
calculation of the ^Be(p,7)^B radiative capture [29j, which 
is important for the standard solar model and neutrino 
physics (see Fig. [8|. We also developed a three-cluster 
extension of the method to describe the Borromean nuclei 
(e.g., ^He and ^^Li). 

2.1.4- Coupled- cluster method 

The coupled-cluster method [30l [SH |32l |33] exhibits a 
favorable scaling that grows polynomially with the mass 
number of the nucleus and the size of the model space. 
The UNEDF collaboration employed an ?n-scheme-based 
coupled-cluster code [34 and an angular-momentum cou- 
pled code [35 . The latter exploits the preservation of an- 
gular momentum and pushed ab initio computation with 
"bare" interactions from chiral effective field theory [36 to 
medium-mass nuclei [37 . Coupled-cluster theory is based 
on a similarity-transformed Hamiltonian and employs a 
nontrivial vacuum such as the Hartree-Fock state. In prac- 
tice, one iteratively solves a large set of nonlinear coupled 
equations. The exploitation of rotational invariance con- 
siderably reduces the number of degrees of freedom but 
comes at the cost of working in a much more complicated 
scheme (involving angular momentum algebra) that poses 
challenges for a scalable and load-balanced implementa- 
tion. 

During UNEDF, several conceptual advances in physics 
and computing were made with the coupled-cluster 
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Figure 8: Ab initio calculated ^Be{p,j)^B S-factor (solid line) com- 
pared with experimental data and the calculation used in the latest 
evaluation (dashed line); see [29] for details. 



method. On the physics side, these include the angular- 
momentum coupled implementation of the coupled-cluster 
method [37 , the use of a Gamow basis for computation 
of weakly bound nuclei ^SJ |39], a practical solution to 
the center-of-mass problem in nuclear structure compu- 
tations |40j, the extension of the method to nuclei with 
up to two nucleons outside a closed subshell [41^, the ap- 
proximation of three-nucleon forces as in-medium correc- 
tion to nucleon-nucleon forces [U, "iS] |44], and the de- 
velopment of theoretically founded extrapolations in finite 
oscillator spaces [45^. On the computational side, seal- 
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Figure 10: Excitation energies of J^ = 2+ states in Ca isotopes. The 
theoretical results (red squares) agree well with data (black circles) 
and predict a soft subshell closure in ^^Ca. 



ing was improved by a work-balancing approach [46l |47] 
based on MPI and OpenMP such that the model-space 
size has increased from ten oscillator shells at the incep- 
tion of UNEDF p^ to 20 oscihator shells at UNEDF's 
completion [44 . Figure [9] shows how adding the use of 
MPI and OpenMP in V2.0 improved the code's scalability 
to thousands of cores, beyond a few hundred cores in Vl.O 
using MPI only, when calculating the small system of ^^Ca 
in 12 oscillator shells. We note that the number of single- 
particle orbitals grows as the third power with the number 
of oscillator shells and that the number of computational 
cycles - in the coupled-cluster method with singles and 
doubles (CCSD) approximation - grows as v?Qn\ (where 
no and n^ are the numbers of occupied and unoccupied 
single-particle states, respectively). Thus, conceptual and 
algorithmic improvements during UNEDF allowed us to 
solve problems that naively required an increase of com- 
putational cycles by about a factor 4,000. The combined 
efforts culminated in the computation of neutron-rich iso- 
topes of oxygen [44 and calcium [49 . 

Doubly magic nuclei are the cornerstones for our under- 
standing of entire regions of the nuclear chart within the 
shell model. For this reason, studies on the evolution of 
structure in neutron-rich semi-magic isotopes of oxygen, 
calcium, nickel, and tin are central to experimental and 
theoretical efforts. With ^^'^^Ca being doubly magic nu- 
clei, many studies were aimed at understanding the struc- 
ture of the rare isotopes ^^'^^Ca and questions regarding 
the AT = 32, 34 sheh closures [50l EH US US |Mj • 

A first-principles description of rare calcium isotopes is 
challenging because it requires the control and understand- 
ing of continuum effects (due to the weak binding) and 
3NFs (as often pivotal contributions arise at next-to-next- 
to leading order in chiral effective field theory [ESI [56l [57] ) . 
Reference [49 reports coupled-cluster results for neutron- 
rich isotopes of calcium that include the effects of the 



continuum and 3NFs (see [58j for complementary calcu- 
lations). It predicts a soft subshell closure in the A^ = 32 
nucleus ^^Ca and an ordering of single-particle orbitals in 
neutron-rich calciums that is at variance with naive shell- 



model expectations. Figure 10 shows the computed ener- 
gies of the first excited J^ = 2+ state in some isotopes of 
calcium and compares them with available data. The high 
excitation energy in ^^Ca is due to its double magicity, 
and the somewhat increased excitation energies in ^^'^^Ca 
suggest that these nuclei exhibit a softer subshell closure. 
Where data are available, the theoretical results agree well 
with experiment. For ^^Ca, theory made a prediction that 
has recently been verified experimentally [59] . 

2.2. Configuration interaction 

The nuclear shell model has been very effective in 
describing the physics of larger nuclei beyond the cur- 
rent reach of pure ab initio methods; indeed, Eugene 
Wigner, Maria Goeppert-Mayer, and J. Hans D. Jensen 
were awarded the 1963 Noble prize for the fundamental 
symmetries and mean field features that underlie the suc- 
cessful nuclear shell model. The shell model for larger 
nuclei uses the same configuration interaction methods as 
the NCSM methods described previously, but with more 
truncated model spaces where not all nucleons are "active" 
and with effective interactions tailored for these spaces. 

Since there are numerous challenging physical appli- 
cations in nuclear physics that vary across the periodic 
table, different CI approaches are needed to efficiently ex- 
ploit the available computational resources. CI approaches 
developed or improved within UNEDF include the follow- 
ing: 

• No-core shell model in the m-scheme basis (mfdu 
[ll[18l|19]; BiGSTiCK [60] El Eg). 



• No-core shell model in a coupled angular momentum 
basis (mfdnj [631164] ). 

• Shell model with a core in a coupled angular momen- 
tum basis (nushellx [65] [62]). 

UNEDF took advantage of common elements in the vari- 
ous CI approaches to improve the effectiveness of the nu- 
clear shell model for all nuclei. 

These CI codes utilize an input NN interaction file and 
a Coulomb interaction between the protons. They all work 
in the neutron-proton basis (i.e., break isospin) and allow 
for charge-dependent NN interactions. In addition, several 
of these codes accept 3NFs as input. All these codes eval- 
uate the spectra, wavefunctions, and a suite of observables 
for low- lying states of the nucleus. 

The implemented algorithms differ considerably among 
the codes as well as support systems for processing the 
output files generated, such as the wavefunctions and one- 
body density matrices, both static and transition. Nu- 
merous cross-comparisons between the codes have been 
accomplished and their respective accuracies confirmed. 
Eigenenergies are obtained to the accuracy of 1 keV or 
better. Other observables are found to differ at the level 
of a few percent because of numerical noise in the wave- 
functions. 

Except for mfdnj (which followed MFDn), the codes 
evolved along independent paths, which emphasized vari- 
ous strategic physics and technological goals. For example, 
the challenges of addressing heavier nuclei impel work- 
ing with a nuclear core; the challenges of working with 
leadership-class machines versus local clusters drive some 
of the algorithmic decisions. The burden of communica- 
tions and memory restrictions help resolve the challenge 
of store-in-memory versus recompute-on-the-fly strategies 
that are implemented differently in these CI codes. 

In light of the need to store large amounts of data for 
retrieval, postanalyses, and reproducibility, we have devel- 
oped a prototype database management system. This pro- 
totype records in the database the metadata of every run. 
The data referenced in the database may include physical 
observables, one-body density matrices, and wavefunctions 
that result from the ab initio codes; such data are typi- 
cally stored on the platforms where runs are performed. 
A user can access this database over the web and find out 
whether the runs of interest have already been performed 
and where the results may be located. 

2. 3. Nuclear density functional theory 

Because of the enormous configuration spaces involved, 
the properties of complex heavy nuclei are best described 
by the superfluid nuclear density functional theory [66 
- rooted in the self-consistent Hartree-Fock-Bogoliubov 
(HFB), or Bogoliubov-de Gennes, problem. The main in- 
gredient of nuclear DFT is the effective interaction be- 
tween nucleons captured by the energy density functional. 
Since the nuclear many-body problem involves two kinds of 



fermions, protons and neutrons, the EDF depends on two 
kinds of densities and currents [671 EH] • isoscalar (neutron- 
plus-proton) and isovector (neutron-minus-proton). The 
coupling constants of the nuclear EDF are usually ad- 
justed to selected experimental data and pseudodata ob- 
tained from ab initio calculations. The self-consistent HFB 
equations allow one to compute the nuclear ground state 
and a set of quasiparticles that are elementary degrees 
of freedom of the system and that can be used to con- 
struct better approximations of the excited states. The 
HFB equations constitute a system of coupled integro- 
differential equations that can be written in a matrix form 
as a self-consistent eigenvalue problem, where the depen- 
dence of the HFB Hamiltonian matrix on the eigenvectors 
(quasiparticle wavefunctions) induces nonlinearities. 

The atomic nucleus is also an open system having 
unbound states at energies above the particle emission 
threshold, and this has implications for the nuclear DFT. 
The finiteness of the HFB potential experienced by a nu- 
cleon implies that the energy spectrum of HFB quasipar- 
ticles contains discrete bound states, resonances, and non- 
resonant continuum states [69 . The size of the continuum 
space may become intractable, especially for complex ge- 
ometries where self-consistent symmetries are broken. To 
this end, one has to develop methods [70^ to treat HFB res- 
onances and nonresonant quasiparticle continuum without 
resorting to the explicit computation of all states. 

The application of high-performance computing, mod- 
ern optimization techniques, and statistical methods has 
revolutionized nuclear DFT during recent years, in terms 
of both developing new functionals and carrying out ad- 
vanced applications. Optimizing the performance of a sin- 
gle HFB run is crucial for making the EDF optimization 
[TTJ [72] manageable and quickly computing tables of nu- 
clear observables [73l [TH [TS] [76] , in order to assess theo- 
retical uncertainties. These advances are described in the 
following sections. 

2.3.1. DFT solvers 

Solutions of HFB equations can be obtained either by 
direct numerical integration on a mesh, provided proper 
boundary conditions are imposed on the domain, or by 
expansion on a basis. For the latter case, the harmonic 
oscillator (HO) basis proves particularly well-adapted to 
nuclear structure problems, as it offers analytical, local- 
ized solutions with convenient symmetry and separability 
features. Although solving the HFB equations for a given 
nuclear configuration is relatively fast on modern comput- 
ers, accurate characterization of nuclear properties often 
requires simultaneous computations of many different con- 
figurations, from a few dozen (e.g., one-quasiparticle con- 
figurations in odd mass nuclei) to a few billion or more in 
extreme applications (such as probing multidimensional 
potential energy surfaces of heavy nuclei during the fission 
process). 

The two primary DFT solvers based on HO expansion 
used by the collaboration are hfbtho [TT and hfodd 



[75] : see [79] and [80], respectively, for their latest re- 
leases. Both codes solve the HFB equations for generalized 
Skyrme functionals in a deformed HO basis and have been 
carefully benchmarked against one another up to the 1 eV 
level. HFBTHO assumes axial and time-reversal symmetry 
of the solutions, making it a very fast program (execu- 
tion completes in typically less than 1 minute on a sin- 
gle node). It is particularly suited for EDF optimization 
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(see Sec. [2.3.3 ) or large-scale surveys of nuclear properties 



[74l[75]. The solver HFODD is fully symmetry-unrestricted: 
this versatility is necessary for science applications such as 
the computation of fission pathways [81 or description of 
high-spin states [82 . 

The new versions of each solver benefited significantly 
from recent advances in high-performance computing and 
from collaborations with computer scientists in UNEDF. 
By expanding the use of tuned BLAS and lapack libraries, 
significant performance gains were reported for both codes 
and enabled new, large-scale studies [83 . The speed of 
HFBTHO was further improved by a factor of 2 by incor- 
porating multithreading; HFODD was turned into a hybrid 
MPI/OpenMP program: nuclear configurations are dis- 
tributed across nodes, while on-node parallelism is imple- 
mented via OpenMP acceleration. 

Figure [TT] illustrates two algorithmic improvements to 
the DFT solver HFODD. The implementation of the Broy- 
den method for nonlinear iterative problems [85 has re- 
duced substantially the number of iterations needed to 
converge the solution in practical applications. The sec- 
ond example shows the application of the augmented La- 
grangian method (ALM) to fission in ^^^Fm [84J. This 
method is generally used for constrained optimization 
problems; it allows precise calculations of multidimen- 
sional energy surfaces in the space of collective coor- 
dinates. Indeed, while the standard quadratic penalty 
method often fails to produce a solution at the required 
values of constrained variables on a rectangular grid, the 
ALM performs well in all cases. Both improvements dis- 
played in Fig. [TT] are key to producing realistic large-scale 
surveys of fission properties in heavy nuclei on leadership- 
class computers, where walltime is limited and expensive. 

Another HFB solver developed by UNEDF is hfb- 
AX. It is based on the B-splines representation of co- 
ordinate space and preserves axial symmetry and space 
inversion [86 . The solver has been carefully bench- 
marked with HFBTHO and used in several applications 
involving complex geometries, such as fission [87 and 
competition between normal superfluidity and Larkin- 
Ovchinnikov (LOFF) phases of polarized Fermi gases in 
extremely elongated traps [88]. Hybrid parallel program- 
ming (MPI+OpenMP) has been implemented in hfb-AX 
to treat large box sizes that are important for weakly 
bound heavy nuclei. 

New generations of DFT solvers will be taking advan- 
tage of emerging architectures, such as GPUs, and new 
programming paradigms. In particular, the cost of per- 
forming dense linear algebra in both hfbtho and HFODD 
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Figure 11: Algorithmic improvements to HFODD. Top: Convergence 
for a typical HFB calculation in the ground state of ^^^Dy with 
HFODD version 2.49t 80 . Using the Broyden method to iterate the 
nonlinear HFB equations has provided significant acceleration com- 
pared with traditional linear mixing techniques. Bottom: Compari- 
son between the augmented Lagrangian method (black squares) and 
the standard quadratic penalty method (open squares) for the con- 
strained HFB calculations of the total energy surface of ^^-^Fm in a 
two-dimensional plane of elongation, Q20, and refiection-asymmetry, 
Qso. (From [84].) 



can become prohibitive as the size of the HO basis in- 
creases, especially for more realistic energy functionals in- 
volving some form of nonlocality; this necessitates novel 
techniques to handle many-body matrix elements [89]. 
The massive amount of data generated by large-scale DFT 
simulations will also require significant investments in vi- 
sualization and data-mining techniques. 

2.3.2. Multiresolution 3D DFT framework 

A parallel, adaptive, pseudospectral-based solver, 
MADNESS-HFB, has been developed to tackle the fully 
symmetry-unrestricted HFB problem for both real and 
complex wave functions in large and asymmetric boxes. 
The main mathematical and algorithmic advantage of 
MADNESS-HFB is its multiscale-multiresolution and sparse 
approximation of functions and the application of oper- 
ators in coordinate space with guaranteed accuracy but 
finite precision, madness-hfb prefers to work with func- 
tions and operators with pseudo-spectral approximations 
based on a multiwavelet basis (up to order 30). Since 
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Figure 12: Quasiparticle wavefunction for a DFT simulation (left, 
top) and its six levels of mult iresolut ion structure (left, bottom). 
The refinement structure is especially noticeable at levels 5 and 6. 
Right: The parallel speedup of one iteration of madness-hfb, for 
solving the DFT problem for 1,640 3D quasiparticle wavefunctions 
with over 4.4 billion equations and unknowns; this simulation was 
performed within a box with a spatial dimension of 120 fermis, using 
8 multiwavelets, up to level 8+ of refinement, and with a relative 
precision of 10~^. 



the multiwavelets consist of smooth, singular, and discon- 
tinuous functions with spatial locality (compact support), 
they are well suited for localized approximation of weak 
singularities and discontinuities or regions of high curva- 
ture ^90j |9i| i92j. Gibbs effects are also reduced. The 
object-oriented (00) nature of the software and template- 
based programming allow each wavefunction and each in- 
tegral or differential operator to have its own boundary 
condition and its own sparse pseudospectral expansion. 
The usual boundary conditions (e.g., Dirichlet, Neumann, 
Robin, quasi-periodic, free, and asymptotic conditions) are 
supported. Fast applications of Green's function for the 
direct solution of Poisson's equation and the Yukawa scat- 
tering kernel are available [93l [SU [95] . In the multiwavelet 
representation, these approximate Green's functions and 
their applications are again based on sparse data with 
guaranteed precision, in contrast to dense tensors based 
on the use of some other basis sets. Other Green's func- 
tions can also be constructed. 

If desired, the user can specify solvers and routines 
from other dense and sparse linear algebra packages such 
as LAPACK or SCALAPACK. For example, parallel and vec- 
torized adaptive quadrature permit the construction of 
the Hamiltonian matrix in the usual manner by using 
the ^2 norm. The Hamiltonian can be diagonalized by 
using multithreaded LAPACK (or a parallel eigensolver), 
and the eigenvectors can be converted back to coefficients 



for the multiwavelet representation. Other capabilities, 
such as high-order approximation of propagators and time- 
stepping required for the solution of time-dependent DFT, 
are also available from applications in time-dependent 
molecular DFT, as well as from simulation of attosecond 
dynamics [96", "97] . 

Underlying this mathematical capability is a parallel 
runtime system that permits the software to scale to hun- 
dreds of thousands of processors and runs on platforms 
from laptops to leadership-class computers. The ability 
to use laptops and workstations is particularly attractive 
for model and code development and testing. In addi- 
tion, the embedding of a parser permits the 00-based 
C++ templated codes representing operations on the co- 
efficients of each wavefunction to be executed as parallel 
tasks. This parser permits out-of-order, distributed mul- 
tithread executions with task- and data-dependency anal- 
ysis. This reduces the stalling of execution units due to 
data dependencies. A user-configured and executed paral- 
lel load-balancing method is also available, as is a parallel 
checkpoint and restart method. 

The 3D MADNESS-HFB has been benchmarked with the 
sphne-based 2D solver hfb-AX [86], 3D hfodd [80], and 
the ID code hfbrad [98 for a variety of problems. Be- 
cause MADNESS-HFB has no limit on the size of the com- 
putational domain, we were able to capture quasiparticle 
wavefunctions with long tails or nonsymmetric potentials 
with steep curvatures and cut-offs to overcome some of the 
limitations of the other solvers. The adaptive structure is 



illustrated in Fig. 12 



The current MADNESS- hfb approach to the HFB prob- 
lem is as follows [99 . Let the coefficients of the wavefunc- 
tions in the tensor product multiwavelet representation be 
the unknowns. The user provides an initial relative pre- 
cision, a set of initial wavefunctions (e.g., in terms of the 
HO basis, splines, etc.), and boundary conditions to start 
the iterative procedure. All the functions, potentials, op- 
erators, and expansion lengths are adapt ively represented 
as needed by the user-defined precision. A generalized ma- 
trix eigenvalue problem is formed by adaptive quadrature. 
The solution eigenvectors are converted to a sparse mul- 
tiwavelet representation for updating the lengths of the 
expansion and the coefficients in the potentials, gradients, 
and other terms before the next iteration and diagonaliza- 
tion. The speed and performance depend on the number 
of coefficients. Usually, the simulations begin with a low 
relative precision, to capture the low-order terms quickly, 
before adaptively increasing the order of approximation 
and the precision for more accurate results. 

2.3.3. EDF optimization 

One of the focus areas of UNEDF was the development 
of an optimization protocol for determining the coupling 
constants of nuclear EDFs. In particular, the collabora- 
tion paid special attention to estimating the errors associ- 
ated with such a procedure and exploring the correlations 
among the coupling constants. The UNEDF optimization 
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protocol was established by focusing on the Skyrme energy 
density. We recall that, in this framework, the energy of 
an even-even nucleus in its ground state is a functional of 
the one-body density matrix and the pairing tensor. The 
Skyrme energy density reads 
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where the isospin index t labels isoscalar {t=0) and isovec- 
tor (t=l) densities, pt is the one-body density matrix, and 
Tt and Jt are derived from pt [67j. In the pairing chan- 
nel, we took a density-dependent pairing energy density 
with mixed surface and volume nature, characterized by 
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the two pairing strengths Vq ^ and Vq for neutrons and 
protons, respectively. The set of coupling constants C^^ , 
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Vq , and Vq^ are the parameters x to be determined. 



The development of fast DFT solvers (see Sec. 2.3.1), 
together with the availability of leadership-class comput- 
ers, permitted us for the first time to set up an opti- 
mization protocol at a fully deformed HFB level. Our 
first parametrization, unedfO, was obtained by consider- 
ing only three types of experimental data: nuclear binding 
energies of both spherical and deformed nuclei, nuclear 
charge radii, and odd-even mass differences in selected nu- 
clei [71]. After recognizing that deformation properties 
needed to be better constrained [lOOj, a fourth data type, 
corresponding to excitation energies of fission isomers in 
the actinides, was added. The resulting parametrization, 
unedfI, gave a significantly better description of fission 
properties £2^, see Fig. flS] (bottom). With the oncoming 
UNEDF2 parametrization, we will expand the optimization 
data set with single-particle level splittings. The new data 
are expected to better constrain the tensor coupling con- 
stants and improve single-particle properties. 

Formally, we solve the optimization problem 
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where d G M"^^ represents the experimental data, w > 
represent weights, and the parameters x to be determined 
are possibly restricted to lie in a domain fl. This prob- 
lem is made difficult because some of the derivatives with 
respect to the parameters x, \/xSi{x)^ may be unavailable 
for some of the theory simulation observables Si . 

Traditional approaches for solving ^ in the absence 
of derivatives typically either estimate these derivatives 
by finite differencing or treat x^ ^is a black-box function of 
X. The former approach can be sensitive to the choice of 
the difference parameter, and care must be taken that the 
expense of the differencing does not grow unnecessarily as 
the number of parameters Ux grows. The latter neglects 
the structure (in the form of the rid residuals) inherent to 

In UNEDF, we instead employed a new optimization 
solver, POUNDerS, that exploits the structure in nonlinear 
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Figure 13: Top: Performance of the POUNDerS algorithm on the min- 
imization of the x^ of Eq. (121 as compared with the standard Nelder- 
Mead method. Bottom: Fission pathway for ^^^Pu along the mass 
quadrupole moment Q20 calculated with SkM, unedfO, and unedfI 
EDFs. The experimental energy of fission isomer (Ejj) and the inner 
(Ea) and outer (Eb) barrier heights are indicated |72| . 



least-squares problems and avoids directly forming compu- 
tationally expensive derivative approximations. POUNDerS 
follows a model-based Newton-like approach, where the 
first- and second-order information is inferred by itera- 
tively forming local interpolation models for each resid- 
ual. Figure 13 (top) shows the efficiency of the solver: 



not only does it converge faster than the standard Nelder- 
Mead algorithm, but it also gives a more accurate solution. 
POUNDerS is available in the open-source Toolkit for Ad- 
vanced Optimization (TAO [101 ). 

2.34. Neutron droplets and DFT 

The properties of homogeneous and inhomogeneous 
neutron matter play a key role in many astrophysical sce- 
narios and in the determination of the symmetry energy 
[ 102, ilQ3[ 1104] . The equation of state of homogeneous neu- 
tron matter has been studied in many earlier investigations 
(see, e.g., [105"). Since neutron matter is not self-bound, 
inhomogeneous neutron matter has been theoretically in- 
vestigated by confining neutrons in external potentials. 
Although neutron drops cannot be realized in experimen- 
tal facilities, they provide a model to study neutron-rich 
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Figure 14: Calculated total energies for neutron droplets in hw = 
5 MeV and 10 MeV harmonic potentials as a function of the neutron 
number A^. The figure shows AFDMC, GFMC, SLy4, and adjusted 
SLy4 results of 109 together with the unedfO and unedfI results. 



isotopes |1Q6[ I1Q7[ 1108] and can bridge ab initio methods 
and DFT. The external potential confining neutrons has 
been chosen to change the geometry and density of the 
system. A Woods-Saxon form produces saturation, mak- 
ing neutron drops similar to ordinary nuclei. Instead, a 
harmonic potential permits one to better control the cal- 
culation of larger systems and to test the approach to the 
thermodynamic limit. 

Nuclear EDFs are commonly optimized to reproduce 
properties of nuclei close to stability, with close numbers 
of protons and neutrons. The use of such functionals to 
study neutron-rich nuclei or the neutron star crust requires 
large extrapolations in neutron excess. In [109] . neutron 
droplets were studied by using QMC methods starting 
from a realistic nuclear Hamiltonian that includes the Ar- 
gonne AV8' two-body interaction supported by the Ur- 
bana IX three-body force. This Hamiltonian fits nucleon- 
nucleon phase shifts, gives a satisfactory description of 
light nuclei, and produces an equation of state of neu- 
tron matter that is compatible with recent neutron star 
observations [110 . The neutron drop's energy calculated 
by using QMC methods was compared with DFT calcu- 
lations. The QMC results showed that commonly used 
Skyrme EDFs typically overbind neutron drops and that 
this effect is due mainly to the neutron density gradient 
term. The adjustment of the gradient together with the 
pairing and spin-orbit terms improves the agreement be- 
tween ab initio QMC calculations with Skyrme both for 
the energy and for neutron densities and radii |109j . 

These results can be compared with the predictions of 
unedfO and unedfI EDFs. Figure [M] shows the calcu- 
lated total energies for neutron droplets in hw = 5 MeV 
and 10 MeV harmonic potentials. The auxiliary field dif- 
fusion Monte Carlo (AFDMC) and GFMC QMC results 
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Figure 15: Top: Deformation energy curves for ^^^Zr calculated 
using microscopic EDFs derived from chiral EFT interactions at dif- 
ferent orders [111) . Bottom: Comparison of microscopic EDF calcu- 
lations of neutron drops at increasing levels of approximation with 
full NCFC calculations starting from the same Hamiltonian 112 . 



of [109 , calculated with the AV8'H-UIX interactions, agree 
well with the DFT calculations [72 . These are encourag- 
ing, since neither unedfO nor unedfI was optimized to 
the pure neutron matter data. Future EDF optimization 
schemes will use ab initio results on neutron droplets as 
pseudodata to improve EDF properties in very neutron- 
rich nuclei. 

2.S.5. Ah initio junctionals 

In parallel with efforts to improve the optimization 
of nuclear EDFs with conventional Skyrme- type terms, 
UNEDF members sought to construct ab initio func- 
tionals based on microscopic chiral effective field theory 
(EFT) |113j . A pathway to such functionals was opened 
with the development of new renormalization group meth- 
ods, which led to softer nuclear Hamiltonians, including 
three-body forces [114 . These soft interactions dramat- 
ically improve convergence properties in many-body cal- 
culations [115 , extending the reach of ab initio methods 
to heavier systems [116^ I117[ I118J . At the same time, 
they make feasible the construction of a microscopically 
based EDF using many-body perturbation theory [119j 
together with improved density matrix expansion (DME) 
techniques |120[ I12H I122i . Carrying out this long-term 
program by individual researchers would be a formidable 
task, but progress was made possible within UNEDF by 
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teaming up with two of the physics-CS/AM partnerships 
described earher. 

An intermediate step toward a fuhy ab initio EDF was 
a new hybrid functional that incorporated long-range chi- 
ral EFT interactions to describe pion-range physics and a 
set of Skyrme-like contact interactions with coupling con- 
stants to be fit. The resulting functional has a much richer 
set of density dependencies than do conventional Skyrme 
functionals. These were incorporated in the DFT solvers, 
and new preoptimization procedures were developed by 
the DFT functional group |lllj . A proof-of-principle test 
in the top panel of Fig. [15] shows deformation energies 
in ^^^Zr calculated using the DME functional at differ- 
ent orders in the chiral expansion (LO, NLO, N^LO). The 
deviations from the Skyrme result show nontrivial effects 
from the finite-range nature of the underlying NN and 3N 
interactions [111 . On-going work includes a rigorous op- 



timization with the procedure outlined in Sec. 2.3.3 and 
then detailed evaluations of the predictive power of the 
DME functional. 

In order to directly validate the new DME procedures 
used in [1111, it was necessary to benchmark against ex- 
act results. The first-ever such calculations were made 
possible by teaming up with the NCSM-MFDn effort (see 
Sec. 2.1.2) using neutron droplets as a controlled the- 
oretical test environment as in Sec. 12.3.41 The DME 
functional was constructed and evaluated for the same 
(model) Hamiltonian used to generate exact results from 
MFDn 11121 for different numbers of neutrons and varied 



traps. Figure [15] (bottom) shows the agreement between 
no-core full configuration (NCFC) results and microscopic 
EDF calculations at different levels of approximation |112j , 
which validates the optimal strategy used to construct 
a microscopically based EDF (the points labeled "fit"), 
while establishing theoretical error bars. Further impor- 
tant DME developments made by external collaborators 
in the FIDIPRO project [123,124 wih be tested in future 
investigations. 

24. Beyond DFT 

Static DFT provides excellent tools for investigating 
nuclear binding energies and other ground-state proper- 
ties. In certain cases, it also can be used to treat dy- 
namical processes. The path to scission during fission, for 
example, sometimes can be predicted accurately by static 
DFT. A reliable description of excitation/decay and reac- 
tions, however, usually requires methods that go beyond 
static DFT. Since an ab initio treatment of the nuclear 
time-evolution is difficult, we employ extensions of DFT 
and related ideas. The simplest extension, the quasiparti- 
cle random phase approximation (QRPA), can be viewed 
as an adiabatic approximation to the linear response in 
time-dependent DFT. It provides the entire spectrum of 
excitations with the same EDF used in static DFT. The 
adiabatic approximation is, of course, severe (as are the 
approximations in the density functional itself) but can 
be applied in any nucleus and folded with reaction theory. 



DFT-based QRPA and its applications to nuclear excita- 
tion and reactions are discussed in Sec. 12.4.11 

DFT-based methods that go beyond the adiabatic ap- 
proximation are also now in use. One can exploit the 
relatively simple dynamics of Fermi gas systems to con- 
struct an approximate time-dependent extension of DFT, 
the time-dependent superfluid local density approxima- 
tion (TDSLDA). The approximation and related compu- 
tational techniques can be applied to such classic problems 
as photoabsorption but also to other time-dependent pro- 
cesses that go beyond linear response. The TDSLDA and 
its applications are discussed in Sec. |2.4.2[ 

We also need efficient methods to accurately compute 
average properties of excited states, such as spin- and 
parity-dependent level densities, which suffice to treat re- 
actions that proceed primarily through a compound nu- 
cleus. Obtaining these densities through a direct diago- 
nalization of the nuclear Hamiltonian and a subsequent 
level counting is not efficient, but several techniques based 
on statistical spectroscopy can be used instead. However, 
even statistical spectroscopy poses computational chal- 
lenges that demand high-performance computational tech- 
niques and resources. Some advances in computational 
spectroscopy, leading to the first accurate calculation of 
densities of levels with unnatural parity, are described in 
Sec. EMI 

2.4' F QRPA and reactions 

Members of the UNEDF collaboration developed and 
exploited both an extremely accurate spherical Skyrme 
QRPA code [125^ and an equally accurate, though com- 
putationally much more intensive, deformed (axially sym- 
metric) Skyrme QRPA code [126'. The latter, which can 
treat both spherical and deformed nuclei, is at the fore- 
front of the modern QRPA. Other groups have developed 
their own versions of the deformed Skyrme, Gogny, or rel- 
ativistic QRPA [127, UMl [HH USUI [m]; most of these 
have some disadvantages compared with ours (e.g., a lack 
of full self-consistency, oscillator bases that don't capture 
continuum physics, etc.) but also the occasional advan- 
tage (e.g., full continuum wavefunctions rather than the 
approximate representation of the continuum we describe 
below) . 

Both our spherical and deformed codes diagonalize 
the traditional QRPA A-B matrix [132 , constructed from 
single-quasiparticle states in the canonical basis [132 in a 
large box (typically 20 fm in each coordinate), so that con- 
tinuum states are taken into account in discretized form. 
Both codes work with arbitrary Skyrme density function- 
als plus delta pairing, include all rearrangement terms, and 
break neither parity nor time-reversal symmetries. Both 
output transition amplitudes to the entire spectrum of ex- 
cited states. 

The two codes have some differences as well. The 
spherical code gets its single-quasiparticle wavefunctions, 
represented on an equidistant mesh, from an HFB program 
called HFBMARIO, which derives from the code HFBRAD 
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[69] . The deformed code takes its wavefunctions from the 
Vanderbilt HFB program |133j , which uses B-sphnes to 
represent wave functions. Each QRPA code represents 
those wavefunctions in the same manner as the HFB pro- 
grams it rehes on. 

Both QRPA codes have been tested in many ways, in- 
cluding against one another. With the spherical code, we 
calculated energy-weighted sums in Ca, Ni, and Sn iso- 
topes from the proton drip line through the neutron drip 
line for J^ = 0+, 1~, and 2+ multipoles with Skyrme 
parameter sets SkM* and SLy4 and found excellent agree- 
ment with analytical values [125 . Spurious states in the 
J^ = 0+, 1+, and 1~ channels are well separated from 
physical states in both codes, though the spherical one 
performs a bit better because it can include all combina- 
tions of HFB two-quasiparticle states in the QRPA basis 
without making the calculation intractable. 

The collaboration used the spherical QRPA to study 
systematics of 2+ states across the table of isotopes and 
for microscopic calculations of reaction rates; they used 
the deformed version for a more limited study of 2+ states 
and giant resonances in rare-earth nuclei [134 . 

The collaboration also used transition densities from 
the spherical QRPA to calculate nucleon-nucleus scatter- 
ing. The transition amplitudes produced by our spherical 
matrix QRPA, when combined with single-particle wave 
functions, yield radial transition densities. These can in 
turn be folded with the interaction between the projectile 
and the nuclear constituents (i.e., the nucleon-nucleon in- 
teraction) to produce transition potentials that excite tar- 
get states. References |135[ 1136) 1137] report the develop- 
ment of a code to fold the densities for all QRPA states be- 
low 30 MeV with a Gaussian-shaped nucleon-nucleon po- 
tential. The result is a microscopic coupled-channels cal- 
culation that successfully produces angular distributions 
and inelastic cross sections for nucleon-induced reactions — 
quantities that can be compared directly with scattering 
data — at scattering energies between 10 and 70 MeV. To 
satisfactorily describe observed absorption, we had to ex- 
plicitly couple also to all one-nucleon pickup channels lead- 
ing to intermediate deuteron formation. Figure 16 illus- 
trates the effect of such couplings on nucleon-induced ab- 
sorption cross sections. The direct connection between the 
calculated cross sections and the nuclear structure ingre- 
dients makes this kind of reaction calculation a good test 
of the structure model. 

The collaboration also took significant steps to develop 
a much more efficient implementation of the QRPA. The 
finite amplitude method [ 140, 14 1 allows one to effectively 
take the derivatives of mean fields that enter the QRPA 
equations numerically, through relatively straightforward 
modifications to the mean-field codes themselves. A simple 
iterative procedure then solves the equations. Our initial 
application, to monopole resonances in the deformed nu- 
cleus ^^^Pu [142 , consumes a small fraction of the time our 
matrix QRPA implementation would use (see 1143(1144] for 
complementary work based on iterative Arnoldi diagonal- 
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Figure 16: Reaction cross section as a function of incident energy 
for p + ^"-"Zr. The results are shown for couphngs to the inelas- 
tic states (dash-dotted line) and to the inelastic and transfer chan- 
nels with nonorthogonality corrections (solid line). The Koning- 
Delaroche 138 optical model calculations are also shown (dashed 
line). (Data from [139].) 



ization). 

2.4-2. Time- dependent DFT for superfluid systems 

The application of DFT to nuclear physics requires two 
nontrivial elements: the ability to describe both superflu- 
idity and time-dependent phenomena. In order to avoid 
the nonlocal character of the DFT extension to super- 
fluid systems, the superfluid local density approximation 
(SLDA) and its time-dependent extension TDSLDA have 
been developed [14^ [Tim [T471 [T48l [T49l [THIl [TCTl [15^ [15^ 



SLDA and TDSLDA have been applied to a large num- 
ber of fermionic systems and phenomena: vortex structure 
in neutron matter and cold atomic systems, generation and 
dynamics of quantized vortices and their crossing and re- 
connection, excitation of the Anderson-Higgs modes, the 
LOFF phase, quantum shock waves and excitation of do- 
main walls, one- and two-nucleon separation energies, gi- 
ant dipole resonance in superfluid triaxial nuclei, and com- 
plex collisions. In Fig. [T7| we illustrate the case of a head- 
on collision of two superfluid fermion clouds, which was 
studied experimentally. Both SLDA and TDSLDA are de- 
rived by using appropriately determined EDFs with QMC 
input for homogeneous systems and validating the predic- 
tions on independent QMC calculations of inhomogeneous 
systems in the well-studied case of a unitary Fermi gas; 
see [147, 148, 150^ for details. The form of the EDF for 
a unitary Fermi gas is largely determined by dimensional 
arguments; translational, rotational symmetry, and parity; 
gauge and Galilean covariance (which specifies the depen- 
dence on current densities); and renormalizability of the 
TDSLDA formalism. 

For nuclear systems we lack ab initio results of the same 
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Figure 17: Three consecutive frames of the head-on cohision of two 
fermion clouds of ?^750 particles in which quantum shock waves and 
domain walls/solitons (topological excitations) are formed |152| . The 
X- and y-directions have an aspect ratio of ~30. 



quality and rely on a more phenomenological approach, 
but with significant microscopic input. The nuclear EDFs 
should satisfy the usual symmetries [68 and the consis- 
tency with the best available ab initio results. 

The numerical implementation of the SLDA and TD- 
SLDA equations leads to hundreds of thousands of cou- 
pled nonlinear 3D time-dependent PDEs, which are solved 
by using the discrete variable representation approach 
[T551 [156] on desktops [Mil [Ml [Hi USD] and— as a re- 
sult of UNEDF collaborations with computer scientists — 
leadership-class supercomputers |15Q[ I15H I152[ 11531 1154] . 
In Fig. [18] we illustrate the first calculation of the photoex- 
citation of a triaxial superfluid nucleus performed within 
TDSLDA (^^^Os) and two other axially deformed nuclei, 
as well as a comparison with the absolute experimental 
data (without any fitting parameters). The determination 
of the ground-state properties of these nuclei and their 
subsequent time-evolution required full diagonalizations of 
Hermitian matrices of sizes up to 5 • 10^ x 5 • 10^ and solu- 
tions of 5 • 10^ coupled time-dependent 3D PDEs. Further 
studies of excitation of medium- and heavy-mass nuclei 
with 7-rays, neutrons, relativist ic heavy ions, and induced 
nuclear fission are the next steps. 

2.4-3. Level densities 

The properties of the excited states of nuclei are key to 
reliably describing reactions and decays. One important 
type of reaction mechanism is the compound nuclear reac- 
tion, which can be described with the statistical model 
of Hauser and Feshbach |157j . The important ingredi- 
ent entering the Hauser-Feshbach theory is the spin- and 
parity-dependent nuclear level density (NLD). Experimen- 
tal information about NLD is limited for stable nuclei and 
not available for radioactive nuclides of interest for nu- 
clear astrophysics. Therefore, a large effort is underway to 
accurately calculate NLD, and an interacting shell model 
approach would be the best model taking into account the 
relevant many body correlations beyond DFT. A direct 
approach by direct CI diagonalization and level counting 
is not feasible because of the exponential increase in CI 
dimensions. We recently proposed [158j an approach to 
calculate shell model spin- and parity-dependent NLD us- 
ing methods of statistical spectroscopy. In addition, we 
showed [159 how one can improve this approach to cal- 
culate the unnatural parity NLD by removing the con- 
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Figure 18: Photoabsorption cross section (sohd black hne) calcu- 
lated within TDSLDA using two Skyrme force parameterizations for 
three deformed open-shell nuclei and the experimental (7, n) cross 
sections (solid purple circles with error bars); see 153 for details. 
With dashed (green), dotted (red), and dot-dashed (blue) lines, we 
display the contribution to the cross section arising from exciting the 
corresponding nucleus along various symmetry axes. 



tribution due to the spurious center-of-mass excitations. 
The associated algorithms were implemented in a high- 
performance computer code, jmoments, [160. 1161} I162J , 
which runs on massively parallel computers and scales well 
up to 10,000 processors [16011162]. 



Figure 19 shows positive- and negative-parity NLD for 
^^Al calculated with jmoments compared with the avail- 
able experimental data obtained by level counting. Some 
known levels have no clear assignment of the parity, which 
leads to upper and lower limits. The calculated positive- 
parity NLD is not new, an sd-shell calculation being avail- 
able for some time. However, the negative-parity NLD was 
calculated only recently by our approach [161 . 

3. Uncertainty Quantification 

Uncertainty quantification is a key element for assess- 
ing the predictive power of a model. When working with 
effective theories with degrees of freedom relevant to the 
problem, the parameters of the theoretical model often 
need to be adjusted to the empirical input. To quantify the 
model uncertainties, sensitivity analysis yields the stan- 
dard deviations and correlations of the model parameters, 
usually encoded as a covariance matrix [165. .166^ .1671 [T68] . 

The calculation of the covariance matrix requires com- 
puting derivatives of the observables with respect to the 
model parameters. When a closed- form expression for the 
derivatives is not available, we estimate the derivatives nu- 
merically using finite differences. To account for the nu- 
merical uncertainty associated with the underlying DFT- 
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Figure 19: Nuclear-level densities for positive parity (red curve) and 
negative parity (black curve) of ^^Al compared with experimental 
data; the solid and dotted staircases represent upper and lower limits, 
respectively. Positive parity NLD is larger than negative parity NLD. 



based calculations, we compute the "noise level" of each 
observable following the approach in ^169] . The difference 
parameters used for estimating the Jacobian matrix asso- 
ciated with ([2| are then obtained using these noise levels 

[TTni . 

Uncertainty quantification was one of the key topics of 
the EDF optimization work performed in the UNEDF col- 



laboration [711 [72]. The upper panel of Fig. 20 shows the 
UNEDF 1 correlation matrix, obtained from the sensitivity 
analysis. As can be seen, some of the surface parame- 
ters of the UNEDF 1 EDF are strongly correlated. In [76 
we used this information to assess the robustness of the 
current EDFs in the predictions of the nuclear landscape 
limits. This is illustrated in the lower panel of Fig. [20j 
which shows calculated and experimental two-neutron sep- 
aration energies for the isotopic chain of even-even zir- 
conium isotopes. The differences between model predic- 
tions are small in the region where data exist and grow 
steadily when extrapolating toward the two-neutron drip 
line (S2n — 0)- Nevertheless, the consistency between the 
models was found to be surprisingly good. This study re- 
quired massive parallel calculations of the nuclear mass 
tables [75]. 

4. High-Performance Computing Resources 

UNEDF science has benefited from access to some of 
the largest computers in the world, provided primarily by 
DOE's Innovative and Novel Computational Impact on 
Theory and Experiment (INCITE) program [171^ . In par- 
ticular, the largest computations of UNEDF were carried 
out on the "Jaguar" machine at Oak Ridge National Lab- 
oratory and the "Intrepid" machine at Argonne National 
Laboratory. Jaguar has gone through several processor 
upgrades during the project, taking it from 30,976 cores 
(Cray XT4 in 2008) to 298,592 cores (Cray XK6 in 2012); 
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Figure 20: Top: UNEDF 1 correlation matrix. Presented are the ab- 
solute values of the correlation coefficients between the parameters 
characterizing the energy density (ITl. Bottom: Theoretical extrapo- 
lations toward drip lines for the two-neutron separation energies Sin 
for the isotopic chain of even-even Zr isotopes using different EDFs 
(SLy4, sv-min, unedfO, unedfI) 76 and frdm [163| and hfb-21 
|164| mass models. Detailed predictions around Sin — are illus- 
trated in the inset. The bars on the SV-min results indicate statistical 
errors due to uncertainty in the coupling constants of the functional. 



Intrepid is an IBM Blue Gene/P with 163,840 processing 
cores. 
Figure 
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shows the UNEDF utilization of these 
computing resources over the years 2008-2013 provided 
through INCITE. The figure highlights the increasing de- 
mand for computing time in low-energy nuclear physics 
research. The combined 2008 INCITE utilization across 
Jaguar and Intrepid was nearly 20 million core-hours and 
by 2012 had increased fourfold. This growth illustrates 
the increasing application of high-performance computing 
in nuclear theory enabled by the physics/computer sci- 
ence/applied mathematics collaborations fostered by UN- 
EDF. 

For the 2013 calendar year, members of the SciDAC-3 
NUCLEI project [172 were granted the sixth largest allo- 
cation of the 61 INCITE projects awarded, with a total al- 
location of 155 million core-hours across three leadership- 
class computing resources. Titan, Mira, and Intrepid. Ti- 
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Figure 21: UNEDF allocation and utilization (in millions of core- 
hours) of leadership-class computing resources from 2008 to 2013. 



tan is a Cray XK7, a hybrid CPU-GPU system with 
299,008 CPU cores and 261,632 GPU streaming multipro- 
cessors, and Mira is an IBM Blue Gene/Q with 786,432 
processing cores. The substantial changes to comput- 
ing systems at both Argonne and Oak Ridge, indicative 
of future trends in high-performance computing, create 
new computational challenges but also new possibilities 
to achieve larger and more accurate calculations. Through 
the close collaborations enabled through UNEDF, and now 
NUCLEI, members are working to continuously scale codes 
to increase physics capabilities and improve performance 
for efficient utilization of these leadership-class resources. 

5. Conclusions 

The examples presented here illustrate the multifaceted 
outcomes of the UNEDF project, both in terms of land- 
mark calculations of nuclear structure and reactions and in 
terms of how nuclear theory is done. The project was very 
productive, as can be assessed by going to the project's 
website, http://unedf.org, which documents the con- 
crete deliverables of UNEDF: publications, highlights, re- 
ports, conference presentations, and computer codes. UN- 
EDF also placed great importance on recruiting the next 
generation of scientists. Annually it provided training to 
30 young researchers. The UNEDF experience has been 
a springboard for advancement, with many UNEDF post- 
docs obtaining permanent positions at U.S. universities, 
national laboratories, and overseas institutions. 

By fostering broad new collaborative efforts between 
physicists, mathematicians, and computer scientists, the 
SciDAC-2 UNEDF project showed how to tackle scientific, 
algorithmic, and computational challenges in the era of 
extreme-scale scientific computing. This effort continues 
with the SciDAC-3 NUCLEI project [172], which builds 
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on the successful strategies of UNEDF. Figure 22 shows 
the key elements of NUCLEI. 
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